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Introduction-Experimental Setup 


T-rays: appropriateness for imaging applications 


Terahertz rays (T-rays) are pulses of sub- picosecond duration of 
electromagnetic radiation. They lie in the region on the electromagnetic 
spectrum between what is traditionally considered electronics and 
photonics. In practical implementations, picosecond bursts may be 
artificially generated in such a way so that they attain near linear phase and 
a broad fractional bandwidth in the Terahertz (THz) frequency range. 
Accordingly, it is possible to directly measure the temporal electric field 
using T- rays in a similar way that this is done for rays that lie in other 
regions on the electromagnetic spectrum. This allows probing depths of 
materials by looking at the arrival time of transmitted or reflected pulses. 
Note that one THz is equal to 1012 Hertz which corresponds to sub 
millimeter wavelengths in free space. Indeed, the wavelengths within the 
bandwidth of a typical T-ray are between 1 mm at 0.3 THz and 0.3 mm at 1 
THz and evidently this can be translated into high temporal resolution (.3 
mm in transmission and .15 in reflection), in imaging and detection 
applications. 
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In addition, T-rays have a number of unique material responses. On the one 
hand, plastics, papers, and many packaging materials are virtually 
transparent to T-rays. In fact, plastics work with T-rays similarly to glass 
with visible light. As an example, all of the lenses that are used in practice 
to suitably focus T-rays are made out of plastic and thus they can be readily 
machined. On the other hand, metals, such as Aluminum, are highly 
reflective. Water is strongly absorbing and T-rays cannot penetrate it, 
rendering biomedical imaging using T-rays of limited interest. However, T- 
rays can be used in package inspection or for manufacturing quality control. 
In 1995, Hu and Nuss took the first T-ray pictures of a semiconductor (see 
references). This work fueled interest in using T-rays for imaging in a 
number of applications and in a variety of configurations. NASA is using T- 
rays to inspect the foam on the shuttle tanks for defects, which is believed 
to have caused the Columbia accident. 


Experimental setup that provided the data used in this project 


The present project uses experimental data obtained by a reflection 
computed tomography method utilizing T-rays (TRCT), proposed by J. 
Pearce, H. Choi, D. Mittleman, J. White, D. Zimdars (see references). In 
TRCT, the objective is to reconstruct the reflectivity edge map of an 
object’s thin tomographic slice. The setup of the TRCT imaging system is 
shown in the next figure. A THz transceiver is used in order to generate T- 
ray pulses and measure the back-reflected waveforms. A lens is placed in 
front of the transceiver to collimate the beam and then the cylindrical lens 
of focal length 12 cm focuses the beam in the vertical direction, confining it 
to a horizontal plane. The generated T-rays illuminate a thin cross-sectional 
Slice of the object of interest. The object is placed on a rotation stage. It is 
rotated in 1 degree increments for 360 degrees and a measurement is taken 
at each viewing angle. 
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The available experimental data for the purposes of this project concern the 
metal square post shown in the next picture. It is a square inch in 
dimensions and made of aluminum. The reason for the choice of this test 
object is that it contains both large and small features. The dents in the 
metal are less than 1 mm across and it is interesting to verify the spatial 
resolution capabilities of the T-rays to resolve these features. Aluminum is a 
strong reflector so the measured reflected waveform is substantial. Since 


aluminum is opaque, it is not expected that the screw hole will be resolved 
because none of the T-rays will reach it. 
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Two main steps for imaging the test object: I) Deconvolution, 
IT) Reconstruction 


The techniques and signal processing algorithms that are used are 
completely analogous to reflection computed tomography used with 
ultrasonic rays (see references). It is evident from the experimental setup 
that reflections from points that lie on the wavefront will sum together. This 
corresponds to a parallel beam projection. Each reflected waveform is the 
convolution of the incident T-ray pulse with the projections of the object’s 
edge map. 


The plan of attack for image retrieval consists of 2 steps. First, the 
projections have to be deconvolved from the measured waveforms. Second, 


the reflectivity map of the object is obtained from the projections. Herein, 
in order to deconvolve the projections from measurements, two methods 
were attempted: Regularized inverse filtering and Wiener filtering. Next, 
the image retrieval procedure is accomplished by application of the Filtered 
Backprojection Algorithm. FBP filters each projection with a ramp filter 
and then backproject the filtered projections across image plane. 


Description/Manipulation of Data 
The raw data at our disposal consisted of: 


e 360 measured waveforms, one from each angle. Each waveform is a 
back- reflection of a transmitted THz pulse at different angles. 

e The reference pulse. This is the received signal reflected off a mirror 
placed at the center of the turning table. We need to know the shape of 
the emitting pulses in order to factor out its contribution from the 
measured waveforms. 

¢ The background signal. This is obtained by just receiving the noise of 
the environment. No pulse is emitted and no specimen is placed on the 
turning table. 


All the above temporal signals were sampled with time step 0.016945 pico- 
seconds (ps) for 400 ps duration. 


In the first plot of figure 1 an example waveform (measured at 70 degrees) 
is plotted against the reference pulse and the background signal. Clearly, 
both the waveform and the pulse follow the trend of the background signal. 
This suggests that the background signal should be first subtracted from the 
reference pulse and all waveforms, before any other signal processing takes 
place. The second plot of the same figure shows the reference and the 
example waveform after subtracting the environmental interference. It 
should be noted that the background signal is mainly a result of random 
mechanical vibration of the transceiver and experiment setup, plus the 
environmental noise that always exists. 


An important observation on the second plot is that the reference pulse 
attains its maximum value just before 200 ps. This shows that it takes 
nearly 200 ps for the THz pulse to travel from the transceiver to the middle 
of the table and back. Moreover, the example waveform reaches the first 
local maximum at around 120 ps, which identifies the closest edge of the 
specimen. 


The initial part of the example waveform consists of only high frequency 
and very low magnitude noise. This is the case for all the measured 
waveforms, although the length of the initial noisy part differs. The noise 


exhibited in the initial parts is representative of the noise present throughout 
all the waveform observations. Thus, it may serve as data to estimate the 
variance of the noise. The estimated variance for all the waveforms 
assuming Gaussian white noise ranges from 10-6~10-7. This verifies that 
the level of noise is very low as can be seen in the example waveform and 
in the reference pulse. 

Reference and example waveforms 


Raw data 
O4 T T T T T T 


Observed measurement at angle 70 degrees 
background signal 
reference pulse al 


Electric field (arb. units) 


200 
Time delay (ps) 


Data after subtracting the background signal 


Observed measurement at angle 70 degrees 
reference pulse | 


Electric field (arb. units) 
S 
nN 


200 
Tire delay (ps) 


Deconvolution with Inverse and Weiner Filters 


Problem Statement 


The first task is to deconvolve the reference pulse from all the measured 
waveforms to obtain the actual projections of our specimen. In general the 
deconvolution problem can be modeled as in the following figure. 

Block diagram of our system 
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In the above figure, p(t) is the projection of the object slice that is to be used 
in the object reconstruction. The reference pulse is denoted by r(t) and 
interpreted as a low- pass degradation filter that convolves with the 
projections. The result of this convolution plus an additive noise n(t) is the 
measured waveforms s(t). In this case, both r(t) and s(t) are known and the 
goal is to obtain p(t) by deconvolving r(t) from s(t) and denoising. 
Mathematically the above block diagram can be written as: 


S(@,f) = [r(t—2u/ c) p,(u)du 


Inverse Filter 


Most of deconvolution schemes are implemented in the frequency domain. 
This makes sense since convolution in the time domain is mapped as 


multiplication in the frequency domain. Inverse filtering is the simplest and 
most naive method for deconvolution. No denoising takes place in this case. 


Equation 1 can be written in the frequency domain as: 


S(@)=Riw)P(@) 


where S(@), R(@) and P(@) are the Fourier transforms of s(t), r(t) and p(t) 
respectively. 
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The main drawback of the method is that R(@) is usually a low-pass filter 
and therefore 1/R(@) is a high-pass filter which attains large values as the 
frequency increases. Thus, equation (3) becomes numerically unstable for 
small R(@) values and greatly amplifies the high- frequency noise 
contribution. This makes inverse filtering very sensitive to even the small 
amounts of high frequency noise which exists in the measured waveforms. 


P(@)= 


One method of decreasing the noise sensitivity inherent when inverse 
filtering is to bound the frequency response 1/R(@) to some threshold y as 
follows: 
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The code that implements inverse filtering is inverseFilter.m 


Wiener Filter 


One of the most widely used restoration techniques is the Wiener filter. 
Contrary to the inverse filtering this method also attempts to diminish noise 
while restoring the original signal. It executes an optimal balance between 
inverse filtering and noise smoothing in the mean square error sense. 
Assuming white Gaussian noise, the Wiener filter expressed in the Fourier 
domain is written as: 


— _R(@)"Spo() 
| R(@) ? Spp(@) + on? 


where Spp(q@) is the power spectrum of the input projection and o2 is the 
variance of the Gaussian noise. The derivation of this formula is based on a 
stochastic framework and is beyond the scope of this project. 


In this case the projections are unknown and an estimate of Spp is needed. 


This is provided by noting that Sss(@)= Spp(@)*|R(@)|2, solving for Spp(@) 
and substituting into equations (5) to get: 
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Interestingly, the last equation can be interpreted as two different filters in 
cascade in the frequency domain. The first is denoising, while the second is 
exactly the inverse filter considered in the previous section. 


The code that implements Wiener filtering is wienerFilter.m . 


Driver for filters 


Results of Deconvolution 


For the purposes of this project, much experimentation was performed in 
order to define the value for the threshold y that gives the best results for 
both inverse and Wiener filtering methods. In addition, much effort was put 
forth to obtain a meaningful estimate for the variance of the noise under the 
Gaussian white noise assumption. To this point, one criterion for the 
effectiveness of each method is subjective observation of the deconvolved 
waveforms as opposed to the measured ones. The first plot in the following 
figure is the previously used example measured waveform. The next two 
plots show the deconvolved waveform (estimate of projection at angle 70), 
using the inverse and Wiener filter algorithms with appropriate parameters. 
Original, inverse filtered, and weiner filtered signals 
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Certain conclusions can be drawn from the above plots. At first, it is evident 
that the patterns of the reference THz pulse signal are attenuated from the 
waveform after the deconvolution process in both cases. This makes the 
main two pulses of the measured waveform signal (occurring at about 120 
and 175 ps), more distinct from the rest of the signal’s features. Intuitively, 
this will result in a sharper image after reconstruction. Secondly, in the case 
of the Wiener filtered signal the absence of high frequency noise is obvious. 
However, since this method balances between inverse filter and noise 
reduction, the inverse filtering part is weaker. Overall, just inverse filtering 


the waveforms turns out to give better results and by appropriately choosing 
the y constant the noise level is kept at negligible level. Note that after 
deconvolution the absolute magnitude of the signal is greatly reduced, but 
this is not an issue for the purpose of reconstruction since relative 
reflectivity values are used. 


A second global criterion is a two dimensional plot of the measured and the 
deconvolved waveforms as a function of delay and angles. These sorts of 
plots are often called sinograms because of the sinusoidal behavior of the 
features. Each sinogram has 360 columns so that each column is a temporal 
plot of the estimated projections. 

Sinograms of original, inverse filtered, and weiner filtered signals 
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Both deconvolved sinograms exhibit sharper main sinusoid patterns, while 
the ripples behind them are less evident than those observed in the 
measured waveform sinogram. 


Reconstruction 


Theory of Filtered Backprojection Algorithm (FBP) 


The FBP algorithm allows us to take the projections, PO(t), developed in 
the previous sections and reconstruct the original image, f(x,y). 


A key idea is the Fourier Slice Theorem. It says that the Fourier Transform 
of a projection at an angle theta is equivalent to the values of the 2- 
dimensional Fourier Transform of the image evaluated along a radial line of 
the same angle. Knowing this fact, we are able to obtain the Fourier 
Transform of the image, F(u,v), from projections taken at multiple angles. 


We start with 2-dimensional Inverse Fourier Transform: 
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Since we have projections for given angles, a change to polar coordinates is 
useful. 
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Using symmetry, this simplifies to: 


fay)= [ F@.Olole”™ dodo 


Using the Fourier Slice Theorem, we substitute in the Fourier Transform of 
the projection, SO(@). 
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With this formula, we are now able to reconstruct the original image. We 
now see that that the FPB algorithm has certain benefits. We can begin 
reconstructing the image after the first projection has been calculated, since 
the image is built up by summing over all the angles. This could increase 
speed and practicality for real time applications. 


Backprojection Implementation 


Our implementation uses three pieces of Matlab code: shrinkPR, filtersinc, 
and backproject. 


ShrinkPR 

The original projection data was greatly over-sampled, taking over 20,000 
samples in 400 picoseconds for each of the 360 angles. In order to reduce 
computation time, we uniformly reduced the number of samples and the 
number of angles. Our shrinkPR code takes the original PR matrix, a matrix 
whose columns are the projections, and outputs a smaller PR matrix with 


fewer samples and angles. It also creates a column vector, theta, which 
contains the angles of the corresponding projections. 


Filtersinc 


This code takes in the PR matrix and outputs the filtered PR matrix. 


O,= | S,(alole?*#do 


—o 


It creates the ramp filter |@| and multiplies it with the FFT of each 
projection. It then takes the IFFT of each filtered projection. 


Backproject 


Backproject implements the discrete approximation of : 
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This code backprojects each of the filtered projections over the image plane 
and sums them together to produce the final reconstructed image. In order 
to implement the continuous function given discrete data points, the 
“round” function was used, effectively interpolating the data. 


Representative Results 
The following are reconstructed images: 


1500 samples per projection and 360 angles: 


reconstructed image 


From this image you can easily see the original object. 


1100 samples per projection and 360 angles: 


Reconstructed Image 


This image is not as clear since fewer samples were used. 


1500 samples and angles 1-180. 


Reconstructed Image 


From this image you can see the left half of the object since we used the 
first half of the projections. 


Conclusions and References 


Conclusions 


Using the terahertz reflected waveforms, we were able to measure the 
projections and reconstruct the original image. This process was completed 
in two steps. In the first, inverse and Wiener filtering were used to 
deconvolve the data from the reference pulse to obtain the actual 
projections of the test object. In the second, the Filtered Backprojection 
Algorithm featuring the Fourier Slice Theorem was used to filter the 
projections using a ramp filter and backproject the result over the image 
plane, thus reconstructing the image of the object. 


As far as the deconvolution part of the project concerns, the regularized 
inverse filter gives better results than Wiener filtering, as already pointed 
out. Care should be exercised when picking the regularized parameter y, so 
as not to increase the noise level of the resulting signal. Usually, this is a 
case-dependent procedure that takes numerical experimentation. It should 
be noted that the original data at hand were of very good quality with low 
noise level. Thus the results of Wiener filtering were worse than those 
obtained by inverse filtering. 


For the reconstruction part, it was found that the number of projection 
angles used was vital to the clarity of the final image. We were able to 
greatly downsample the data and still maintain image quality to a certain 
point. Due to the size of the data, the algorithms ran for minutes. 


Future Work 


Much potential for future improvement and implementation is possible 
using this method of computerized tomography. Advanced deconvolution 
methods featuring wavelets for efficient noise reduction can be used for 
isolating the projections out of the measured waveforms. Furthermore, it 
seems reasonable to cut-off the first part of each measured waveform since 
it only contains noise and no useful information for the image 
reconstruction. This can be accomplished by appropriately windowing the 
raw measurements before any other manipulations takes place. Due to the 


linear nature of the process, different algorithm code could have been 
written to start reconstructing the image immediately after the first 
projection had been calculated. This and other efficiency tools could be 
implemented to greatly increase the speed of the overall reconstruction, 
making the process applicable for real-time use. 
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